MODULE hydro_procedures

!Version: 02/07/2008 (version 2.7)

!-------------------------------------------------------------------------------
! Section 1.0   Module Header
!-------------------------------------------------------------------------------
!
! DESCRIPTION: 
!   This module contains all the subroutines used by DigiBog Hydro.
!   
! CURRENT CODE OWNERS:
!   Paul J. Morris, Andy Baird*, Lisa R. Belyea and Ruth Boogert
!   *School of Geography, 
!   University of Leeds, 
!   Leeds
!   LS2 9JT
!   pmorris@mcmaster.ca
!   a.j.baird@leeds.ac.uk
!   l.belyea@qmul.ac.uk
!   r.boogert@qmul.ac.uk
!         

!  ***** BEGIN LICENSE BLOCK *****
!  Version: MPL 1.1
! 
!  The contents of this file are subject to the Mozilla Public License Version
!  1.1 (the "License"); you may not use this file except in compliance with
!  the License. You may obtain a copy of the License at
!  http://www.mozilla.org/MPL/
! 
!  Software distributed under the License is distributed on an "AS IS" basis,
!  WITHOUT WARRANTY OF ANY KIND, either express or implied. See the License
!  for the specific language governing rights and limitations under the
!  License.
! 
!  The Original Code is DigiBog Hydro v2.7.
! 
!  The Initial Developer of the Original Code is
!  Andy Baird & Paul J Morris.
!  Portions created by the Initial Developer are Copyright (C) 2008
!  the Initial Developer. All Rights Reserved.
! 
!  Contributor(s): ----
! 
!  ***** END LICENSE BLOCK *****

! MODIFICATION HISTORY: 
!   Programmer           Date           Modifications
!   ============================================================================
!   Andrew J. Baird      08/04/2005     Original 1.5-D code (v. 1.3)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       24/02/2006     Conversion to 2.5 D and DigiBog Fortran 
!                                       Standards implemented (v. 2.0)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       27/02/2006     Testing, minor corrections (v.2.1)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       05/03/2006     Subroutines 'column_activation' and 
!                                       'steady_state_check' written (v. 2.2)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       19/06/2006     Testing completed (v. 2.3)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       24/04/2007     Further cleaning, including removal of
!                                       replicated spatial step reference in
!                                       'move_water' subroutine (v 2.5)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       09/05/2007     Above-ground storage facilitated in 
!                                       2.5-D version (v 2.6)
!   ----------------------------------------------------------------------------
!   Paul J. Morris       02/07/2008     Final code cleaning, full annotation.
!                                       (v. 2.7)
!   ----------------------------------------------------------------------------
!
! CODE DESCRIPTION: 
!   Language:           Fortran 95 
!   Software Standards: "UK Meteorological Office Fortran 90 Standards".
!   See:
!   http://research.metoffice.gov.uk/research/nwp/numerical/fortran90/
!          f90_standards.html

! DECLARATIONS

  IMPLICIT NONE

! Global type definition:
! Define a real kind type q with at least 8 decimal digits and an exponent range
! from 10**30 to 10**(-30):
  INTEGER, PARAMETER :: q = SELECTED_REAL_KIND(P = 8, R = 30)


  CONTAINS


!-------------------------------------------------------------------------------
! Section 2.0   Determine locations and water-table levels of boundaries
!-------------------------------------------------------------------------------

  SUBROUTINE column_activation(x_extent, y_extent, water_table,                &
                               activation_status, active_columns)

! This subroutine determines which columns are boundary columns, which are within 
! the boundary and which are external to it.
!
! MODIFICATION HISTORY: 
!   Programmer           Date           Modifications
!   ============================================================================
!   Paul J. Morris       05/03/2006     'column_activation' and 
!                                       'steady_state_check' subroutines written
!   ----------------------------------------------------------------------------   
!
! DECLARATIONS

  IMPLICIT NONE 

  !Subroutine arguments 
  !Scalar arguments with INTENT(IN): 
  INTEGER, INTENT(IN) :: x_extent, y_extent

  !Scalar arguments with INTENT(INOUT):
  INTEGER, INTENT(INOUT) :: active_columns
 
  !Array arguments with INTENT(INOUT): 
  REAL(KIND=q), DIMENSION(:,:), INTENT(INOUT) :: water_table

  !Array arguments with INTENT(OUT):
  CHARACTER(8), DIMENSION(:,:), INTENT(OUT) :: activation_status

  !Local scalars:
  INTEGER :: x, y
         
!-- END OF SUBROUTINE HEADER ---------------------------------------------------

    !CALCULATIONS

    !Entire activation_status array initialised to 'off':
    activation_status = 'off'

    !Active columns counter initialised to zero
    active_columns = 0

    !Edge of model set to outside boundary:
    DO x = 1, x_extent
      water_table(x, 1)        = -999.0
      water_table(x, y_extent) = -999.0
    END DO

    DO y = 2, (y_extent - 1)
      water_table(1, y)        = -999.0
      water_table(x_extent, y) = -999.0
    END DO



    !All columns with water table of -999.0, including those at the edge of the 
    !model (set above) and any others as read from file "in_wattable.txt", are 
    !switched off.  Any columns bordering these become boundaries, and their 
    !water tables held constant throughout the model run.  All other columns are
    !switched on:

    DO x = 2, (x_extent - 1)

      DO y = 2, (y_extent - 1)

        IF (ABS(water_table(x, y) + 999.0) < 0.0001) THEN

          activation_status(x, y) = 'off'

        ELSE IF (ABS(water_table(x - 1, y) + 999.0) < 0.0001                   &
            .OR. ABS(water_table(x + 1, y) + 999.0) < 0.0001                   &
            .OR. ABS(water_table(x, y - 1) + 999.0) < 0.0001                   &
            .OR. ABS(water_table(x, y + 1) + 999.0) < 0.0001)THEN

          activation_status(x, y) = 'boundary'

        ELSE 

          activation_status(x, y) = 'on'

          active_columns = active_columns + 1

        END IF

      END DO

    END DO

  END SUBROUTINE column_activation



!-------------------------------------------------------------------------------
! Section 3.0   Initial Transmissivity Calculations
!-------------------------------------------------------------------------------

  SUBROUTINE trans_height(x_extent, y_extent, no_layers, activation_status,    &
                         layer_attributes, transmissivity)

! This subroutine calculates the transmissivity for layer’s height above the 
! impermeable base.


! DECLARATIONS

  IMPLICIT NONE 

  !Subroutine arguments 
  !Scalar arguments with INTENT(IN): 
  INTEGER, INTENT(IN) :: x_extent, y_extent
 
  !Array arguments with INTENT(IN): 
  REAL(KIND=q), DIMENSION(:,:,:,:), INTENT(IN) :: layer_attributes

  INTEGER, DIMENSION(:,:), INTENT(IN) :: no_layers

  CHARACTER(8), DIMENSION(:,:), INTENT(IN) :: activation_status

  !Array arguments with INTENT(INOUT):
  REAL(KIND=q), DIMENSION(:,:,:,:), INTENT(INOUT) :: transmissivity

  !Local scalars:
  INTEGER :: x, y, z 
         
!-- END OF SUBROUTINE HEADER ---------------------------------------------------

    !CALCULATIONS
    DO x = 1, x_extent

      DO y = 1, y_extent

        !Ignore columns external to boundary:
        IF(activation_status(x, y)     == 'on'                                 &
          .OR. activation_status(x, y) == 'boundary') THEN
    
          transmissivity(x, y ,1, 1) = layer_attributes(x, y, 1, 1)
        
          transmissivity(x, y, 1, 2) = layer_attributes(x, y, 1, 1)            &
                                     * layer_attributes(x, y, 1, 2)        
    
          DO z = 2, no_layers(x, y)
      
            transmissivity(x, y, z, 1) = transmissivity(x, y, z-1, 1)          &
                                       + layer_attributes(x, y, z, 1)
          
            transmissivity(x, y, z, 2) = transmissivity(x, y, z-1, 2)          & 
                                       + (layer_attributes(x, y, z, 1)         &
                                       * layer_attributes(x, y, z, 2))
     
          END DO

        END IF

      END DO

    END DO

  END SUBROUTINE trans_height



!-------------------------------------------------------------------------------
! Section 4.0   Water Capacity Initialisation
!-------------------------------------------------------------------------------

  SUBROUTINE layer_water_depth (x_extent, y_extent, no_layers,                 &
                                activation_status, layer_attributes,           &
                                layer_storage)

! This subroutine calculates the amount of water (expressed as a depth) capable 
! of being stored in each layer in each column.

! DECLARATIONS

  IMPLICIT NONE 

  !Subroutine arguments  
  !Scalar arguments with INTENT(IN): 
  INTEGER, INTENT(IN) :: x_extent, y_extent

  !Array arguments with INTENT(IN): 
  INTEGER, DIMENSION(:,:), INTENT(IN) :: no_layers

  CHARACTER(8), DIMENSION(:,:), INTENT(IN) :: activation_status
  
  REAL(KIND=q), DIMENSION(:,:,:,:), INTENT(IN) :: layer_attributes
 
  !Array arguments with INTENT(INOUT): 
  REAL(KIND=q), DIMENSION(:,:,:,:), INTENT(INOUT) :: layer_storage

  !Local scalars:
  INTEGER :: x, y, z
   
!-- END OF SUBROUTINE HEADER ---------------------------------------------------


    !Calculations
    DO x = 1, x_extent

      DO y = 1, y_extent

        !Ignore columns external to boundary:
        IF (activation_status(x, y) == 'on'                                    &
           .OR. activation_status(x, y) == 'boundary') THEN
    
        layer_storage(x, y, 1, 1) = layer_attributes(x, y, 1, 1)

        layer_storage(x, y, 1, 2) = layer_attributes(x, y, 1, 1)               &
                                  * layer_attributes(x, y, 1, 3)        

          DO z = 2, no_layers(x, y)

          layer_storage(x, y, z, 1) = layer_storage(x, y, z-1, 1)              & 
                                    + layer_attributes(x, y, z ,1)
         
          layer_storage(x, y, z, 2) = layer_attributes(x, y, z, 1)             &
                                    * layer_attributes(x, y, z, 3)
     
          END DO

        END IF

      END DO

    END DO

  END SUBROUTINE layer_water_depth


!-------------------------------------------------------------------------------
! Section 5.0   Depth-Averaged Hydraulic Conductivity Calculations
!-------------------------------------------------------------------------------

  SUBROUTINE wat_k_mean(x_extent, y_extent, no_layers, activation_status,      &
                        layer_attributes, transmissivity, wk_mean, water_table)
                        
! This subroutine calculates the depth-averaged K below the water table for each 
! column.
    
! DECLARATIONS

  IMPLICIT NONE 

  !Subroutine arguments   
  !Scalar arguments with INTENT(IN):
  INTEGER, INTENT(IN) :: x_extent, y_extent 
         
  !Array arguments with INTENT(IN): 
  REAL(KIND=q), DIMENSION(:,:), INTENT(IN) :: water_table

  REAL(KIND=q), DIMENSION(:,:,:,:), INTENT(IN) :: layer_attributes,            &
                                                  transmissivity
  
  CHARACTER(8), DIMENSION(:,:), INTENT(IN) :: activation_status
  
  INTEGER, DIMENSION(:,:), INTENT(IN) :: no_layers
  
  !Array arguments with INTENT(INOUT): 
  REAL(KIND=q), DIMENSION(:,:), INTENT(INOUT) :: wk_mean

  !Local scalars: 
  INTEGER :: x, y, z

!-- END OF SUBROUTINE HEADER ---------------------------------------------------


    !CALCULATIONS
    DO x = 1, x_extent

      DO y = 1, y_extent
        
        !Ignore columns external to boundary:
        IF (activation_status(x, y)     == 'on'                                &
           .OR. activation_status(x, y) == 'boundary') THEN
      
          !Water table at/above the peatland surface
          IF (water_table(x, y)                                                &
              >= transmissivity(x, y, no_layers(x, y), 1)) THEN 
        
            wk_mean(x, y) = transmissivity(x, y, no_layers(x, y), 2)           &
                          / transmissivity(x, y, no_layers(x, y), 1)
        

          !Water table below the top of the first layer:
          ELSE IF (water_table(x, y) < transmissivity(x, y, 1, 1)) THEN 
      
            wk_mean(x, y) = layer_attributes(x, y, 1, 2)
        

          !Water table equal to or above top of first layer and below top 
          !of mire surface:
          ELSE  
      
            DO z = 1, (no_layers(x, y) - 1)
            
              !If water table equal to height of layer z:
              IF (water_table(x, y) == transmissivity(x, y, z, 1)) THEN 
          
                wk_mean(x, y) = transmissivity(x, y, z, 2)                     &
                              / transmissivity(x, y ,z ,1)
          

              !If water table above top of layer z:
              ELSE IF (water_table(x, y) > transmissivity(x, y, z ,1)) THEN 
            
                !If water table below top of layer z+1:
                IF (water_table(x, y) < transmissivity(x, y, z+1, 1)) THEN 

                  wk_mean(x, y) = (transmissivity(x, y, z, 2)                  &
                                + (layer_attributes(x, y, z+1, 2)              &
                                * (water_table(x, y)                           &
                                - transmissivity(x, y, z, 1))))                &
                                / water_table(x, y)                      
               
                END IF
            
              END IF
                     
            END DO
         
          END IF

        END IF

      END DO
          
    END DO        

  END SUBROUTINE wat_k_mean


!-------------------------------------------------------------------------------
! Section 6.0   2-Dimensional Flux Calculations
!-------------------------------------------------------------------------------

  SUBROUTINE move_water(x_extent, y_extent, activation_status, base_altitude,  &
                        spatial_step, water_table, wk_mean, timestep,          &
                        water_change, rainfall)

! This subroutine calculates the net movement of water (expressed as a depth) 
! between columns.  It is assumed a harmonic mean operates between the mean K of
! each column.  The output is a change in depth of water (volume per unit area) 
! to be added to or taken from the previous water-table elevation.
    
! DECLARATIONS

  IMPLICIT NONE   
                        
! Subroutine arguments 
! Scalar arguments with INTENT(IN):
  REAL(KIND=q), INTENT(IN) :: spatial_step, timestep, rainfall
                              
  INTEGER, INTENT(IN) :: x_extent, y_extent
  
! Array arguments with INTENT(IN): 
  REAL(KIND=q), DIMENSION(:,:), INTENT(IN) :: base_altitude, water_table,   &
                                                 wk_mean

  CHARACTER(8), DIMENSION(:,:), INTENT(IN) :: activation_status

! Array arguments with INTENT(INOUT):                                               
  REAL(KIND=q), DIMENSION(:,:), INTENT(INOUT) :: water_change

! Local arrays: 
  REAL(KIND=q), DIMENSION(x_extent, y_extent) :: x_flux, y_flux

! Local scalars:
  INTEGER :: x, y   
                         
!-- END OF SUBROUTINE HEADER ---------------------------------------------------


!   CALCULATIONS

    DO x = 1, x_extent

      DO y = 1, y_extent

!       Ignore columns external to boundary:
        IF (activation_status(x, y) == 'on'                                    &
         .OR. activation_status(x, y) == 'boundary') THEN

!         Volume flux (in x direction) between columns using Dupuit-Forchheimer
!         approximation Harmonic mean is used for mean K between columns:
          x_flux(x, y) = (2 * wk_mean(x, y) * wk_mean(x + 1, y))               &
                       / (wk_mean(x, y) + wk_mean(x + 1, y))                   &
                       * ((water_table(x, y) + water_table(x + 1, y))/2)       &
                       * (base_altitude(x, y) + water_table(x, y)              &
                       - base_altitude(x + 1, y)  - water_table(x + 1, y))     &
                       * timestep 
        END IF
                 
      END DO

    END DO	



    DO y = 1, y_extent

      DO x = 1, x_extent

!       Ignore columns external to boundary:
        IF (activation_status(x, y) == 'on'                                    &
         .OR. activation_status(x, y) == 'boundary') THEN
        
!         Volume flux (in y direction) between columns using Dupuit-Forchheimer
!         approximation Harmonic mean is used for mean K between columns:
          y_flux(x, y) = (2 * wk_mean(x, y) * wk_mean(x, y + 1))               &
                       / (wk_mean(x, y) + wk_mean(x, y + 1))                   &
                       * ((water_table(x, y) + water_table(x, y + 1))/2)       &                        
                       * (base_altitude(x, y) + water_table(x, y)              &
                       - base_altitude(x, y + 1) - water_table(x, y + 1))      &
                       * timestep

        END IF
             
      END DO

    END DO



!     Convert volume into depth of water (in absence of soil matrix):
      DO x = 1, x_extent
        DO y = 1, y_extent
          IF (activation_status(x, y) == 'on') THEN

            water_change(x, y) = (x_flux(x - 1, y) - x_flux(x ,y)              &
                               +  y_flux(x, y - 1) - y_flux(x ,y))             &
                               / (spatial_step ** 2)

            water_change(x, y) = water_change(x, y) + (rainfall * timestep)

          END IF
      
        END DO

      END DO

  END SUBROUTINE move_water


!-------------------------------------------------------------------------------
! Section 7.0   Water-Table Update
!-------------------------------------------------------------------------------

  SUBROUTINE water_table_update(x_extent, y_extent, no_layers,                 &
                                activation_status, layer_attributes,           &
                                layer_storage, water_table, water_change)


! This subroutine updates the position of the water table based on the storage 
! available in layers above or below the water table.
    
! DECLARATIONS

  IMPLICIT NONE  

! Subroutine arguments 
! Scalar arguments with INTENT(IN):
  INTEGER, INTENT(IN) :: x_extent, y_extent 
         
! Array arguments with INTENT(IN): 
  REAL(KIND=q), DIMENSION(:,:,:,:), INTENT(IN) :: layer_attributes,            &
                                                  layer_storage
                                                  
  INTEGER, DIMENSION(:,:), INTENT(IN) :: no_layers

  CHARACTER(8), DIMENSION(:,:), INTENT(IN) :: activation_status
 
! Array arguments with INTENT(INOUT):
  REAL(KIND=q), DIMENSION(:,:), INTENT(INOUT) :: water_change, water_table
         
! Local scalars: 
  REAL(KIND=q) :: test_depth
  INTEGER :: x, y, z, marker
  
!-- END OF SUBROUTINE HEADER ---------------------------------------------------
                                 
!   CALCULATIONS
    DO x = 1, x_extent

      DO y = 1, y_extent

!       Only update water table in columns that are switched on:
        IF (activation_status(x, y) == 'on') THEN
           
          !Locate layer in which water table resides:
          
          IF (water_table(x, y) <= layer_storage(x, y, 1, 1)) THEN

            marker = 1
        
          ELSE
                
            DO z = 1, (no_layers(x, y) -1)

              IF (water_table(x, y) > layer_storage(x, y, z, 1)) THEN

!               If water table is in layer z + 1:
                IF (water_table(x, y) <= layer_storage(x, y, z + 1, 1)) THEN

                  marker = z + 1
                
                  EXIT
            
                END IF
          
              END IF
        
            END DO
        
          END IF
      
      
          !Does water table rise?
          IF (water_change(x, y) > 0.0) THEN 
        
            !Ignore any water-table updates if water table already at surface:
            IF (water_table(x, y) < layer_storage(x, y, no_layers(x, y), 1))   &
                                                                            THEN
                                               
              z = marker
                
              IF (water_table(x, y) == layer_storage(x, y, z, 1)) THEN

                !Water table at top of layer z
                !No storage available - advance to next layer:

                z = z + 1
              
              END IF

              DO
              
                test_depth = layer_storage(x, y, z, 2)                         &
                           * (layer_storage(x, y, z, 1) - water_table(x, y))   &
                           / layer_attributes(x, y, z, 1)
                         
                            
                IF (water_change(x, y) <= test_depth) THEN

                
                  !Water table moves up within layer z:
                  water_table(x, y) = water_table(x, y)                        &
                                    + (water_change(x, y)                      &
                                    / layer_attributes(x, y, z, 3))

                  EXIT
              
                ELSE !Water table moves up to next layer or rises to the surface
              
                  !Has surface been reached?
                  IF (z == no_layers(x, y)) THEN

                    water_table(x, y) = layer_storage(x, y, z, 1)
                    EXIT

                  ELSE
                
                    water_table(x, y) = layer_storage(x, y, z, 1)

                    water_change(x, y) = water_change(x, y) - test_depth

                    z = z + 1
                  
                  END IF

                END IF              

              END DO
                   
            END IF

          END IF


          !Does the water table fall?
          IF (water_change(x, y) < 0.0) THEN !Water table falls
        
            !Ignore any water-table updates if water table already at or 
            !below base of column:
            IF (water_table(x, y) > 0.0 ) THEN

              z = marker
                
              DO 
  
                IF (z == 1) THEN

                  test_depth = layer_storage(x, y, 1, 2) * water_table(x ,y)   &
                             / layer_attributes(x, y, 1, 1)
              
                ELSE
              
                  test_depth = layer_storage(x, y, z, 2) *                     & 
                             (water_table(x, y) - layer_storage(x, y, z-1, 1)) &
                             / layer_attributes(x, y, z, 1)
                           
                END IF
                         
                            
                IF (ABS(water_change(x, y)) <= test_depth) THEN 
                                                              
                  !Water table falls within layer z or to top of layer z-1:
                  water_table(x, y) = water_table(x, y) + ( water_change(x, y) &
                                    / layer_attributes(x, y, z, 3) )

                  EXIT
              
                ELSE !Water table falls to top of layer z-1 (if it exists)
                     !and needs to fall further, provided a lower layer
                     !exists:
             
                    !Has base of column been reached?
                    IF (z == 1) THEN

                      water_table(x, y) = 0.0

                      EXIT

                    ELSE
                
                      water_table(x, y) = layer_storage(x, y, z-1, 1)

                      water_change(x, y) = water_change(x, y) + test_depth

                      z = z - 1
                  
                    END IF

                END IF              

              END DO
          
            END IF

          END IF

        END IF

      END DO

    END DO


  END SUBROUTINE water_table_update
  

!-------------------------------------------------------------------------------
! Section 8.0   Check for steady-state hydrological behaviour
!-------------------------------------------------------------------------------

  SUBROUTINE steady_state_check(x_extent, y_extent, activation_status,         &
                                hydro_status, water_change, active_columns,    &
                                steady_threshold, steady_columns)


! This subroutine checks whether the hydrologcal model has reached steady-state
! behaviour in a critical proportion of columns, and if so notifies the main 
! program.

! DECLARATIONS

  IMPLICIT NONE  

! Subroutine arguments 
! Scalar arguments with INTENT(IN):
  INTEGER, INTENT(IN) :: x_extent, y_extent, active_columns

  REAL(KIND=q), INTENT(IN) :: steady_threshold, steady_columns
  
! Array arguments with INTENT(IN): 
  REAL(KIND=q), DIMENSION(:,:), INTENT(IN) :: water_change

  CHARACTER(8), DIMENSION(:,:), INTENT(IN) :: activation_status
 
! Scalar arguments with INTENT(OUT):
  CHARACTER(9), INTENT(OUT) :: hydro_status
         
! Local scalars: 
  REAL(KIND=q) :: check_value
  INTEGER :: x, y
  
!-- END OF SUBROUTINE HEADER ---------------------------------------------------

!   CALCULATIONS
    check_value = 0.0
    
    DO x = 1, x_extent

      DO y = 1, y_extent

        IF (activation_status(x, y) == "on") THEN

          IF (ABS(water_change(x, y)) <= steady_threshold) THEN

            check_value = check_value + 1.0

          END IF

        END IF

      END DO

    END DO

    IF (check_value >= (steady_columns * active_columns)) THEN

      hydro_status = "steady"

    ELSE

      hydro_status = "transient"

    END IF

  END SUBROUTINE steady_state_check
!-------------------------------------------------------------------------------

  
END MODULE hydro_procedures

  
!-------------------------------------------------------------------------------
! End of Module
!-------------------------------------------------------------------------------
